Load libraries
# https://github.com/pbs-assess/gfvelocities/blob/main/R/make_prediction_grid.R
library(tidyverse)
#> ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
#> ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
#> ✔ tibble 3.1.8 ✔ dplyr 1.0.10
#> ✔ tidyr 1.2.0 ✔ stringr 1.5.0
#> ✔ readr 2.1.1 ✔ forcats 0.5.1
#> Warning: package 'tidyr' was built under R version 4.0.5
#> ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
#> ✖ dplyr::filter() masks stats::filter()
#> ✖ dplyr::lag() masks stats::lag()
library(tidylog)
#>
#> Attaching package: 'tidylog'
#>
#> The following objects are masked from 'package:dplyr':
#>
#> add_count, add_tally, anti_join, count, distinct, distinct_all,
#> distinct_at, distinct_if, filter, filter_all, filter_at, filter_if,
#> full_join, group_by, group_by_all, group_by_at, group_by_if,
#> inner_join, left_join, mutate, mutate_all, mutate_at, mutate_if,
#> relocate, rename, rename_all, rename_at, rename_if, rename_with,
#> right_join, sample_frac, sample_n, select, select_all, select_at,
#> select_if, semi_join, slice, slice_head, slice_max, slice_min,
#> slice_sample, slice_tail, summarise, summarise_all, summarise_at,
#> summarise_if, summarize, summarize_all, summarize_at, summarize_if,
#> tally, top_frac, top_n, transmute, transmute_all, transmute_at,
#> transmute_if, ungroup
#>
#> The following objects are masked from 'package:tidyr':
#>
#> drop_na, fill, gather, pivot_longer, pivot_wider, replace_na,
#> spread, uncount
#>
#> The following object is masked from 'package:stats':
#>
#> filter
library(sp)
library(raster)
#>
#> Attaching package: 'raster'
#>
#> The following object is masked from 'package:tidylog':
#>
#> select
#>
#> The following object is masked from 'package:dplyr':
#>
#> select
#>
#> The following object is masked from 'package:tidyr':
#>
#> extract
library(terra)
#> terra version 1.3.22
#>
#> Attaching package: 'terra'
#>
#> The following object is masked from 'package:dplyr':
#>
#> src
library(devtools)
#> Loading required package: usethis
# Source code for map plots
source_url("https://raw.githubusercontent.com/maxlindmark/marine-litter/main/R/functions/map-plot.R")
#> SHA-1 hash of file is 7dd050eb6d584c4f200cdf9c1d93ae83b8ecdecb
#> Warning: package 'sf' was built under R version 4.0.5
#> Linking to GEOS 3.9.1, GDAL 3.4.0, PROJ 8.1.1; sf_use_s2() is TRUE
#> Spherical geometry (s2) switched off
West pred grid
# Read west coast data
litterw <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/marine-litter/main/data/west_coast_litter.csv")
#> Rows: 371 Columns: 23
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (4): haul_id, ansse_area, trend_area, var
#> dbl (18): plast_a, sup_a, fishery_a, tot_a, plast_b, quarter, st_no, haul_no...
#> lgl (1): balse_area
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
x <- litterw$X
y <- litterw$Y
z <- chull(x, y)
coords <- cbind(x[z], y[z])
coords <- rbind(coords, coords[1, ])
plot(coords[, 1] ~ coords[, 2]) # plot data
sp_poly <- sp::SpatialPolygons(
list(sp::Polygons(list(sp::Polygon(coords)), ID = 1))
)
sp_poly_df <- sp::SpatialPolygonsDataFrame(sp_poly,
data = data.frame(ID = 1)
)
class(sp_poly_df)
#> [1] "SpatialPolygonsDataFrame"
#> attr(,"package")
#> [1] "sp"
class(sp_poly)
#> [1] "SpatialPolygons"
#> attr(,"package")
#> [1] "sp"
plot(sp_poly)
plot(sp_poly_df)
cell_width <- 2 # 2*2 km grid cell
pred_grid <- expand.grid(
X = seq(min(litterw$X), max(litterw$X), cell_width),
Y = seq(min(litterw$Y), max(litterw$Y), cell_width),
year = unique(litterw$year)
)
ggplot(pred_grid %>% filter(year == 2019), aes(X, Y)) +
geom_point(size = 0.1) +
theme_void() +
coord_sf()
#> filter: removed 48,240 rows (83%), 9,648 rows remaining
sp::coordinates(pred_grid) <- c("X", "Y")
inside <- !is.na(sp::over(pred_grid, as(sp_poly_df, "SpatialPolygons")))
pred_grid <- pred_grid[inside, ]
pred_grid <- as.data.frame(pred_grid)
plot_map_west +
geom_point(data = pred_grid, aes(X*1000, Y*1000), size = 0.001, alpha = 0.5) +
facet_wrap(~year, ncol = 3) +
geom_sf(size = 0.1) +
NULL
# Add lat and lon
xy <- as.matrix(pred_grid %>% dplyr::select(X, Y) %>% mutate(X = X*1000, Y = Y*1000))
#> mutate: changed 25,362 values (100%) of 'X' (0 new NA)
#> changed 25,362 values (100%) of 'Y' (0 new NA)
v <- vect(xy, crs="+proj=utm +zone=33 +datum=WGS84 +units=m")
y <- project(v, "+proj=longlat +datum=WGS84")
lonlat <- geom(y)[, c("x", "y")]
pred_grid$lon <- lonlat[, 1]
pred_grid$lat <- lonlat[, 2]
ggplot(filter(pred_grid, year == 2017), aes(lon, lat)) + geom_point()
#> filter: removed 21,135 rows (83%), 4,227 rows remaining
# Add ocean area
# https://stackoverflow.com/questions/34272309/extract-shapefile-value-to-point-with-r
# https://gis.ices.dk/sf/
shape <- shapefile("data/assessment_areas_marine_waters/ANSSE_assessmentareas_20181116.shp")
#> Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS
#> = dumpSRS, : Discarded datum European_Terrestrial_Reference_System_1989 in
#> CRS definition: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000
#> +ellps=GRS80 +units=m +no_defs
plot(shape)
shape2 <- spTransform(shape, crs("+proj=longlat +datum=WGS84 +no_defs"))
plot(shape2)
pts <- SpatialPoints(cbind(pred_grid$lon, pred_grid$lat),
proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs"))
pred_grid$area <- over(pts, shape2)$MarUnitID
ggplot(filter(pred_grid, year == 2017), aes(lon, lat, color = area)) + geom_point()
#> filter: removed 21,135 rows (83%), 4,227 rows remaining
# Add depth
dep_raster <- terra::rast("data/Mean depth natural colour (with land).nc")
crs(dep_raster, proj = TRUE)
#> [1] "+proj=longlat +datum=WGS84 +no_defs"
plot(dep_raster)
pred_grid$depth <- terra::extract(dep_raster, pred_grid %>% dplyr::select(lon, lat))$elevation
ggplot(pred_grid, aes(lon, lat, color = depth*-1)) +
geom_point()
pred_grid$depth <- pred_grid$depth*-1
pred_grid <- pred_grid %>% drop_na(depth)
#> drop_na: removed 492 rows (2%), 24,870 rows remaining
plot_map_west +
theme_plot(base_size = 14) +
geom_point(data = pred_grid, aes(X*1000, Y*1000, color = depth, alpha = area)) +
geom_sf(size = 0.1)
#> Warning: Using alpha for a discrete variable is not advised.
# Save
write_csv(pred_grid, "data/pred_grid_west.csv")
East pred grid
# Read east coast data
littere <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/marine-litter/main/data/east_coast_litter.csv")
#> Rows: 368 Columns: 23
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (4): haul_id, balse_area, trend_area, var
#> dbl (18): plast_a, sup_a, fishery_a, tot_a, plast_b, quarter, st_no, haul_no...
#> lgl (1): ansse_area
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
x <- littere$X
y <- littere$Y
z <- chull(x, y)
coords <- cbind(x[z], y[z])
coords <- rbind(coords, coords[1, ])
plot(coords[, 1] ~ coords[, 2]) # plot data
sp_poly <- sp::SpatialPolygons(
list(sp::Polygons(list(sp::Polygon(coords)), ID = 1))
)
sp_poly_df <- sp::SpatialPolygonsDataFrame(sp_poly,
data = data.frame(ID = 1)
)
class(sp_poly_df)
#> [1] "SpatialPolygonsDataFrame"
#> attr(,"package")
#> [1] "sp"
class(sp_poly)
#> [1] "SpatialPolygons"
#> attr(,"package")
#> [1] "sp"
plot(sp_poly)
plot(sp_poly_df)
cell_width <- 2 # 2*2 km grid cell
pred_grid <- expand.grid(
X = seq(min(littere$X), max(littere$X), cell_width),
Y = seq(min(littere$Y), max(littere$Y), cell_width),
year = unique(littere$year)
)
# ggplot(pred_grid %>% filter(year == 2019), aes(X, Y)) +
# geom_point(size = 0.1) +
# theme_void() +
# coord_sf()
sp::coordinates(pred_grid) <- c("X", "Y")
inside <- !is.na(sp::over(pred_grid, as(sp_poly_df, "SpatialPolygons")))
pred_grid <- pred_grid[inside, ]
pred_grid <- as.data.frame(pred_grid)
# plot_map_east +
# geom_point(data = pred_grid, aes(X*1000, Y*1000), size = 0.001, alpha = 0.5) +
# facet_wrap(~year, ncol = 3) +
# geom_sf(size = 0.1) +
# NULL
# Add lat and lon
xy <- as.matrix(pred_grid %>% dplyr::select(X, Y) %>% mutate(X = X*1000, Y = Y*1000))
#> mutate: changed 93,480 values (100%) of 'X' (0 new NA)
#> changed 93,480 values (100%) of 'Y' (0 new NA)
v <- vect(xy, crs="+proj=utm +zone=33 +datum=WGS84 +units=m")
y <- project(v, "+proj=longlat +datum=WGS84")
lonlat <- geom(y)[, c("x", "y")]
pred_grid$lon <- lonlat[, 1]
pred_grid$lat <- lonlat[, 2]
ggplot(filter(pred_grid, year == 2017), aes(lon, lat)) + geom_point()
#> filter: removed 77,900 rows (83%), 15,580 rows remaining
# Add ocean area
# https://stackoverflow.com/questions/34272309/extract-shapefile-value-to-point-with-r
# https://gis.ices.dk/sf/
shape <- shapefile("data/assessment_areas_marine_waters/BALSE_assessmentareas_20181116.shp")
#> Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS
#> = dumpSRS, : Discarded datum European_Terrestrial_Reference_System_1989 in
#> CRS definition: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000
#> +ellps=GRS80 +units=m +no_defs
plot(shape)
shape2 <- spTransform(shape, crs("+proj=longlat +datum=WGS84 +no_defs"))
plot(shape2)
pts <- SpatialPoints(cbind(pred_grid$lon, pred_grid$lat),
proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs"))
pred_grid$area <- over(pts, shape2)$MarUnitID
ggplot(filter(pred_grid, year == 2017), aes(lon, lat, color = area)) + geom_point()
#> filter: removed 77,900 rows (83%), 15,580 rows remaining
# Add depth
dep_raster <- terra::rast("data/Mean depth natural colour (with land).nc")
crs(dep_raster, proj = TRUE)
#> [1] "+proj=longlat +datum=WGS84 +no_defs"
plot(dep_raster)
pred_grid$depth <- terra::extract(dep_raster, pred_grid %>% dplyr::select(lon, lat))$elevation
ggplot(pred_grid, aes(lon, lat, color = depth*-1)) +
geom_point()
pred_grid$depth <- pred_grid$depth*-1
pred_grid <- pred_grid %>% drop_na(depth)
#> drop_na: removed 12,600 rows (13%), 80,880 rows remaining
plot_map_east +
theme_plot(base_size = 14) +
geom_point(data = pred_grid, aes(X*1000, Y*1000, color = depth, alpha = area)) +
geom_sf(size = 0.1)
#> Warning: Using alpha for a discrete variable is not advised.
# Save
write_csv(pred_grid, "data/pred_grid_east.csv")